Introduction au Reservoir Computing avec ReservoirPy¶

Nathan Trouvain
Inria - Mnemosyne




R4 - 9 novembre 2021

 Sommaire¶

  • Concepts et fonctionnalités clés
  • Chapitre 1 : application à des données simples
  • Chapitre 2 : déployer les capacités génératives
  • Chapitre 3 : l'apprentissage en ligne
  • Chapitre 4 : démonstration d'application - chute de robot
  • Chapitre 5 : démonstration d'application - chants de canaris
  • Aller plus loin : comprendre les hyperparamètres

Concepts et fonctionnalités clés ¶

  • Numpy, Scipy, et seulement Numpy et Scipy
  • Exécution efficace (parallélisation des calculs, optimisation des règles d'apprentissage)
  • Plusieurs méthodes d'apprentissage, online, offline
  • Outils d'aide à l'optimisation des paramètres
  • Documentation: https://reservoirpy.readthedocs.io/en/latest/
  • GitHub: https://github.com/reservoirpy/reservoirpy

Informations générales¶

  • Toutes les données sont des tableaux NumPy
  • Le temps est toujours représenté par la première dimension des tableaux/vecteurs.

 Chapitre 1 : Application du Reservoir Computing a des données simples ¶

Série temporelle de Mackey-Glass

L'équation de Mackey-Glass décrit l'évolution de différent phénomènes biologiques, comme la quantité relative de globules rouges dans le temps.

Cette équation est une équation différentielle retardée :

$$ \frac{dP(t)}{dt} = \frac{a P(t - \tau)}{1 + P(t - \tau)^n} - bP(t) $$

où $a = 0.2$, $b = 0.1$, $n = 10$. Le retard $\tau$ est de $17$ pas de temps. $\tau$ a une forte influence sur le caractère chaotique de la série temporelle modelée. $17$ donne des résultats plutôt chaotiques.

In [3]:
from reservoirpy.datasets import mackey_glass
from reservoirpy.observables import nrmse, rsquare

timesteps = 2510
tau = 17
X = mackey_glass(timesteps, tau=tau)

# rescale between -1 and 1
X = 2 * (X - X.min()) / (X.max() - X.min()) - 1
In [5]:
plot_mackey_glass(X, 500, tau)
  • Pas totalement imprédictible... (non aléatoire)
  • ...mais pas facilement prédictible non-plus (non périodique)
  • Comme la variation du rythme cardiaque, la météo, la bourse...

1.1. Tâche 1: prédire la valeur de la série à un horizon de 10 pas de temps¶

Prédire $P(t + 10)$ connaissant $P(t)$.

Préparation les données¶

In [7]:
from reservoirpy.datasets import to_forecasting

x, y = to_forecasting(X, forecast=10)
X_train1, y_train1 = x[:2000], y[:2000]
X_test1, y_test1 = x[2000:], y[2000:]

plot_train_test(X_train1, y_train1, X_test1, y_test1)

Préparation de l'ESN¶

In [8]:
units = 100
leak_rate = 0.3
spectral_radius = 1.25
input_scaling = 1.0
connectivity = 0.1
input_connectivity = 0.2
regularization = 1e-8
seed = 1234

In [10]:
from reservoirpy.nodes import Reservoir, Ridge

reservoir = Reservoir(units, input_scaling=input_scaling, sr=spectral_radius, 
                      lr=leak_rate, rc_connectivity=connectivity, 
                      input_connectivity=input_connectivity, seed=seed)

readout   = Ridge(1, ridge=regularization)

esn = reservoir >> readout
In [11]:
y = esn(X[0])  # initialisation
reservoir.Win is not None, reservoir.W is not None, readout.Wout is not None
Out[11]:
(True, True, True)
In [12]:
np.all(readout.Wout == 0.0)
Out[12]:
True

Entraînement de l'ESN¶

L'apprentissage est offline ("hors-ligne") : il n'a lieu qu'une seule fois, sur l'ensemble des données d'entraînement.

In [13]:
esn = esn.fit(X_train1, y_train1)

2000it [00:00, 12340.04it/s]
In [15]:
plot_readout(readout)

Test de l'ESN¶

In [16]:
y_pred1 = esn.run(X_test1)

100%|███████████████████████████████████████| 500/500 [00:00<00:00, 9522.81it/s]
In [17]:
plot_results(y_pred1, y_test1)

Coefficient de détermination $R^2$ et erreur quadratique normalisée :

In [18]:
rsquare(y_test1, y_pred1), nrmse(y_test1, y_pred1)
Out[18]:
(0.9998810506313958, 0.0026034595249913566)

 1.2 Compliquons la tâche¶

Passons d'un horizon de prédiction de 10 pas de temps à un horizon de 100 pas de temps

In [19]:
x, y = to_forecasting(X, forecast=100)
X_train2, y_train2 = x[:2000], y[:2000]
X_test2, y_test2 = x[2000:], y[2000:]

plot_train_test(X_train2, y_train2, X_test2, y_test2)
In [21]:
plot_results(y_pred2, y_test2, sample=400)

Determination coefficient $R^2$ and NRMSE:

In [22]:
rsquare(y_test2, y_pred2), nrmse(y_test2, y_pred2)
Out[22]:
(0.94803604102033, 0.05593776748240906)

 Chapitre 2 : Déployer les capacités génératives des ESNs ¶

  • Entraînement de l'ESN sur une tâche de prédiction de 1 pas de temps.
  • Test de l'ESN sur ses propres prédictions (mode génératif).
In [23]:
units = 500
leak_rate = 0.3         # - leaking rate
spectral_radius = 0.99  # - rayon spectral
input_scaling = 1.0     # - facteur de mise à l'échelle des entrées (input scaling)
connectivity = 0.1      # - densité des connexions du reservoir vers lui même
input_connectivity = 0.2  # et des entrées vers le reservoir
regularization = 1e-4   # - coefficient de régularisation (L2)
seed = 1234             # reproductibilité

 Entraînement à la prévision sur un horizon court¶

In [25]:
esn = reset_esn()

x, y = to_forecasting(X, forecast=1)
X_train3, y_train3 = x[:2000], y[:2000]
X_test3, y_test3 = x[2000:], y[2000:]

esn = esn.fit(X_train3, y_train3)

2000it [00:00, 8100.12it/s]

 Génération¶

  • 100 pas de temps de la série de test utilisés comme "échauffement";
  • 300 pas de temps générés "à partir de rien";
In [26]:
seed_timesteps = 100

warming_inputs = X_test3[:seed_timesteps]

warming_out = esn.run(warming_inputs, reset=True)  # échauffement

100%|███████████████████████████████████████| 100/100 [00:00<00:00, 4247.18it/s]
In [27]:
nb_generations = 400

X_gen = np.zeros((nb_generations, 1))
y = warming_out[-1]
for t in range(nb_generations):  # génération
    y = esn(y)
    X_gen[t, :] = y
In [28]:
X_t = X_test3[seed_timesteps: nb_generations+seed_timesteps]
plot_generation(X_gen, X_t, nb_generations, warming_out=warming_out, 
                warming_inputs=warming_inputs, seed_timesteps=seed_timesteps)

Chapitre 3 : Apprentissage hors-ligne ¶

Apprentissage se déroulant de manière incrémentale.

Utilisation de l'algorithme FORCE (Sussillo and Abott, 2009)

In [30]:
from reservoirpy.nodes import FORCE

reservoir = Reservoir(units, input_scaling=input_scaling, sr=spectral_radius, 
                      lr=leak_rate, rc_connectivity=connectivity, 
                      input_connectivity=input_connectivity, seed=seed)

readout   = FORCE(1)


esn_online = reservoir >> readout

Entraînement pas à pas¶

In [31]:
outputs_pre = np.zeros(X_train1.shape)
for t, (x, y) in enumerate(zip(X_train1, y_train1)): # pour chaque pas de temps de la série :
    outputs_pre[t, :] = esn_online.train(x, y)
In [32]:
plot_results(outputs_pre, y_train1, sample=100)
In [33]:
plot_results(outputs_pre, y_train1, sample=500)

Entraînement sur une séquence complète¶

In [35]:
esn_online.train(X_train1, y_train1)

pred_online = esn_online.run(X_test1)  # Wout est maintenant figée

100%|███████████████████████████████████████| 500/500 [00:00<00:00, 8039.78it/s]
In [36]:
plot_results(pred_online, y_test1, sample=500)

Determination coefficient $R^2$ and NRMSE:

In [37]:
rsquare(y_test1, pred_online), nrmse(y_test1, pred_online)
Out[37]:
(0.9998104511137069, 0.0032864753423817064)

 Chapitre 4 : Application réelle : chute de robot ¶

In [39]:
features = ['com_x', 'com_y', 'com_z', 'trunk_pitch', 'trunk_roll', 'left_x', 'left_y', 
            'right_x', 'right_y', 'left_ankle_pitch', 'left_ankle_roll', 'left_hip_pitch',
            'left_hip_roll', 'left_hip_yaw', 'left_knee', 'right_ankle_pitch',
            'right_ankle_roll', 'right_hip_pitch', 'right_hip_roll',
            'right_hip_yaw', 'right_knee']

prediction = ['fallen']
force = ['force_orientation', 'force_magnitude']
In [44]:
plot_robot(Y, Y_train, F)

 Entraînement de l'ESN¶

In [46]:
from reservoirpy.nodes import ESN

reservoir = Reservoir(300, lr=0.5, sr=0.99, input_bias=False)
readout   = Ridge(1, ridge=1e-3)
esn = ESN(reservoir=reservoir, readout=readout, workers=-1)  # version distribuée
In [47]:
esn = esn.fit(X_train, y_train)

100%|███████████████████████████████████████| 1444/1444 [00:51<00:00, 27.86it/s]
In [48]:
res = esn.run(X_test)

100%|█████████████████████████████████████████| 361/361 [00:14<00:00, 25.08it/s]
In [51]:
plot_robot_results(y_test, res)
In [52]:
print("RMSE moyenne :", f"{np.mean(scores):.4f}", "±", f"{np.std(scores):.5f}")
print("RMSE moyenne (avec seuil) :", f"{np.mean(filt_scores):.4f}", "±", f"{np.std(filt_scores):.5f}")

RMSE moyenne : 0.1693 ± 0.10356
RMSE moyenne (avec seuil) : 0.1445 ± 0.15254

 Chapitre 5 : Application réelle : le chant des canaris domestiques ¶

Les données peuvent être téléchargées sur Zenodo : https://zenodo.org/record/4736597

In [53]:
im = plt.imread("./static/canary.png")
plt.figure(figsize=(5, 5)); plt.imshow(im); plt.axis('off'); plt.show()
In [55]:
display(audio)
Your browser does not support the audio element.

Plusieurs motifs temporels répétitifs differents à décoder : les phrases.

  • Un label par type de phrase à classifier dans le temps.
  • Un label SIL annotant le silence, à détecter également pour permettre la segmentation du chant.
In [56]:
im = plt.imread("./static/canary_outputs.png")
plt.figure(figsize=(15, 15)); plt.imshow(im); plt.axis('off'); plt.show()

 Entraînement de l'ESN¶

In [61]:
esn = esn.fit(X_train, y_train)

100%|███████████████████████████████████████████| 90/90 [00:09<00:00,  9.22it/s]
In [62]:
outputs = esn.run(X_test)

100%|████████████████████████████████████████| 10/10 [00:00<00:00, 29289.83it/s]
In [64]:
scores  # pour chaque chant testé
Out[64]:
[0.9506604506604507,
 0.9087872559095581,
 0.9531759739013625,
 0.9421525697165245,
 0.9662828947368421,
 0.9536231884057971,
 0.8867638129933212,
 0.9411598302687412,
 0.9209350292196631,
 0.9014840543100726]
In [65]:
print("Précision moyenne :", f"{np.mean(scores):.4f}", "±", f"{np.std(scores):.5f}")

Précision moyenne : 0.9325 ± 0.02503

 Pour aller plus loin¶

  • Essayer de varier le nombre de chants entraînés : combien faut-il de chants au minimum pour reproduire les résultats obtenus sur 100 chants ?
  • Trouver le meilleur jeu d'hyperparamètres : outils basés sur la bibliothèque hyperopt (https://github.com/hyperopt/hyperopt).
  • Essayer d'appliquer un input scaling différent pour chaque groupe de variables d'intérêt.

Merci de votre attention.¶

Nathan Trouvain
Inria - Mnemosyne




R4 - 9 novembre 2021

Aller plus loin : Comprendre les hyperparamètres et leurs effets ¶

In [66]:
units = 100             # - nombre de neurones dans le reservoir
leak_rate = 0.3         # - leaking rate
spectral_radius = 1.25  # - rayon spectral
input_scaling = 1.0     # - facteur de mise à l'échelle des entrées
connectivity = 0.1      # - densité des connexions du reservoir vers lui même
input_connectivity = 0.2  # et des entrées vers le reservoir
regularization = 1e-8   # - coefficient de régularisation (L2)
seed = 1234             # reproductibilité

1. Le rayon spectral¶

Le rayon spectral est la valeur propre maximale de la matrice des poids du réservoir ($W$).

In [67]:
states = []
radii = [0.1, 1.25, 10.0]
for sr in radii:
    reservoir = Reservoir(units, sr=sr, input_scaling=0.1, lr=leak_rate, rc_connectivity=connectivity,
                         input_connectivity=input_connectivity)
   
    s = reservoir.run(X_test1[:500])
    states.append(s)
In [68]:
units_nb = 20

plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
    plt.subplot(len(radii)*100+10+i+1)
    plt.plot(s[:, :units_nb], alpha=0.6)
    plt.ylabel(f"$sr={radii[i]}$")
plt.xlabel(f"Activations ({units_nb} neurons)")
plt.show()

  • $-$ rayon spectral $\rightarrow$ dynamiques stables

  • $+$ rayon spectral $\rightarrow$ dynamiques chaotiques

Rayon spectral et Echo State Property : rayon spectral $\rightarrow$ 1 (assure que les états internes ne sont pas affectés par l'initialisation).

2. Le facteur de mise à l'échelle des entrées (input scaling)¶

Il s'agit d'un coefficient appliqué à $W_{in}$, venant changer l'échelle des données en entrée.

In [69]:
states = []
scalings = [0.01, 0.1, 1.]
for iss in scalings:
    reservoir = Reservoir(units, sr=spectral_radius, input_scaling=iss, lr=leak_rate,
                          rc_connectivity=connectivity, input_connectivity=input_connectivity)
   
    s = reservoir.run(X_test1[:500])
    states.append(s)
In [71]:
units_nb = 20

plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
    plt.subplot(len(scalings)*100+10+i+1)
    plt.plot(s[:, :units_nb], alpha=0.6)
    plt.ylabel(f"$iss={scalings[i]}$")
plt.xlabel(f"Activations ({units_nb} neurons)")
plt.show()

Correlation moyenne des activités des neurones du reservoir avec les entrées :

In [72]:
for i, s in enumerate(states):
    corr = correlation(states[i], X_test1[:500])
    print(f"ISS : {scalings[i]}, correlation moyenne : {corr}")

ISS : 0.01, correlation moyenne : -0.2752749355697493
ISS : 0.1, correlation moyenne : -3.07462701783721
ISS : 1.0, correlation moyenne : -11.893465621922127
  • $+$ input scaling $\rightarrow$ activités corrélées aux données
  • $-$ input scaling $\rightarrow$ activités libres

L'input scaling peut aussi être utilisé pour ajuster l'influence de chaque donnée en entrée.

3. Le taux de fuite de la mémoire (leaking rate)¶

$$ x(t+1) = \underbrace{\color{red}{(1 - \alpha)} x(t)}_{\text{état actuel}} + \underbrace{\color{red}\alpha f(u(t+1), x(t))}_{\text{données suivantes}} $$

avec $\alpha \in [0, 1]$ et:

$$ f(u, x) = \tanh(W_{in} \cdotp u + W \cdotp x) $$
In [73]:
states = []
rates = [0.02, 0.2, 0.9]
for lr in rates:
    reservoir = Reservoir(units, sr=spectral_radius, input_scaling=input_scaling, lr=lr,
                          rc_connectivity=connectivity, input_connectivity=input_connectivity)
   
    s = reservoir.run(X_test1[:500])
    states.append(s)
In [74]:
units_nb = 20

plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
    plt.subplot(len(rates)*100+10+i+1)
    plt.plot(s[:, :units_nb] + 2*i)
    plt.ylabel(f"$lr={rates[i]}$")
plt.xlabel(f"States ({units_nb} neurons)")
plt.show()
  • $+$ leaking rate $\rightarrow$ faible inertie, faible mémorisation des états précédents
  • $-$ leaking rate $\rightarrow$ forte inertie, grande mémorisation des états précédents

Le leaking rate contrôle la "mémoire" de l'ESN. Il peut être vu comme l'inverse de sa constante de temps